#ifndef DREW_INPUT_H_
#define DREW_INPUT_H_
/*
 * Input.h
 *
 *  Created on: Mar 27, 2013
 *      Author: drewlewis
 */

#include <chemistry/molecule/molecule.h>
#include <chemistry/qc/basis/integral.h>
#include <chemistry/qc/basis/split.h>
#include <chemistry/qc/lcao/soad.h>
#include "Global.h"
#include "iostream"

// Returns how big each tile should be along one direciton.
template<typename T, typename S>
std::vector<T>
get_tile_size(sc::Ref<sc::GaussianBasisSet> &basis, std::vector<S> &order){

    std::vector<T> tilesize;
    int trange = 0;
    tilesize.push_back(0);
    //Loop over clusters order[0] = 0
    for(std::size_t i = 1; i < order.size(); ++i){

        //Loop over atoms in cluster
        for(std::size_t j = order[i - 1]; j < order[i]; ++j){
            //For each atom in cluster sum up number of basis functions
            trange += basis->nbasis_on_center(j);
        }

        //Add the number of basis functions to the previous number and pushback
        //This gives functions per cluster
        tilesize.push_back(tilesize[i - 1] + trange);
        trange = 0;
    }

    return tilesize;
}

// Returns a vector of pairs where each pair has the fisrt and last basis
// Funciton on a cluster.
template<typename T, typename S>
std::vector<std::pair<T, T> >
get_cluster_range(std::vector<T> &index, S &clusters){

    std::vector<std::pair<T, T> > range;

    // loop over clusters and return a pair for each cluster with the first and
    // last function on that cluster.
    for(std::size_t i = 0; i < clusters; ++i){
        range.push_back( std::make_pair(index[i], index[i + 1]) );
    }

    return range;
}

//Returns an array that tells tiled array how big to make each tile
inline std::array<TRange1, 2>
rank2_tiled_range(const std::vector<unsigned int> &size){

    typedef TiledArray::TiledRange1 range1;

    std::array<range1, 2> blocking = {{range1(size.begin(), size.end()),
                                       range1(size.begin(), size.end()) }};
    return blocking;
}

//Tells tiled array the blocking information for rank 3 tiles.
inline std::array<TRange1, 3>
rank3_tiled_range(const std::vector<unsigned int> &size1,
                  const std::vector<unsigned int> &size2){

    typedef TiledArray::TiledRange1 range1;
    std::array<range1, 3> blocking = {{ range1(size1.begin(), size1.end()) ,
                                        range1(size1.begin(), size1.end()) ,
                                        range1(size2.begin(), size2.end()) }};
    return blocking;
}

//Tells tiled array the blocking information for rank 4 tiles.
inline std::array<TRange1, 4>
rank4_tiled_range(const std::vector<unsigned int> &size){

    typedef TiledArray::TiledRange1 range1;
    std::array<range1, 4> blocking = {{ range1(size.begin(), size.end()) ,
                                        range1(size.begin(), size.end()) ,
                                        range1(size.begin(), size.end()) ,
                                        range1(size.begin(), size.end()) }};
    return blocking;
}

// Get the molecule information from the input file
inline sc::Ref<sc::Molecule>
get_input_mol(int argc, char** argv, sc::Ref<sc::KeyVal> keyval){

    // get Molecule
    sc::Ref<sc::Molecule> mol_in;
    mol_in << keyval->describedclassvalue("molecule");
    assert(mol_in->point_group()->char_table().order() == 1);

    return mol_in;
}

// Get a regular basis set given a molecule
inline sc::Ref<sc::GaussianBasisSet>
get_obs_basis(sc::Ref<sc::Molecule> &mol){
    sc::Ref<sc::GaussianBasisSet> obs;
    sc::Ref<sc::AssignedKeyVal> akv = new sc::AssignedKeyVal;
    akv->assign("molecule", mol.pointer());
    akv->assign("name", "cc-pVDZ");
    obs = new sc::GaussianBasisSet(sc::Ref<sc::KeyVal>(akv));
    akv->assign("basis", obs.pointer());
    return obs;
}


// Get a density fitting basis set given a molecule
inline sc::Ref<sc::GaussianBasisSet>
get_dfbs_basis(sc::Ref<sc::Molecule> &mol){
    sc::Ref<sc::GaussianBasisSet> dfbs;
    sc::Ref<sc::AssignedKeyVal> akv = new sc::AssignedKeyVal;
    akv->assign("molecule", mol.pointer());
    akv->assign("name", "cc-pVDZ-RI");
    dfbs = new sc::GaussianBasisSet(sc::Ref<sc::KeyVal>(akv));
    akv->assign("basis", dfbs.pointer());
    return dfbs;
}

//Get the soad density guess from mpqc
inline sc::Ref<sc::SymmSCMatrix>
get_soad_density(sc::Ref<sc::Molecule> &mol, Basis &basis){
    sc::Ref<sc::AssignedKeyVal> akv = new sc::AssignedKeyVal;
    akv->assign("molecule", mol.pointer());
    akv->assign("basis", basis.pointer());
    sc::Ref<sc::SuperpositionOfAtomicDensities> soad_guess =
             new sc::SuperpositionOfAtomicDensities(sc::Ref<sc::KeyVal>(akv));
    sc::Ref<sc::SymmSCMatrix> P_guess = soad_guess->ao_density();
    return P_guess;
}

//Read clustering information from input.
std::size_t
inline get_nclusters(std::string filename){

    std::ifstream input;
    input.open(filename.c_str());
    std::size_t nclusters;
    input >> nclusters;
    input.close();
    std::cout << "nclusters = " << nclusters << std::endl;
    return nclusters;
}

inline sc::Ref<sc::Integral>
get_integral_factory(int argc, char** argv){
    sc::Ref<sc::Integral> integral = sc::Integral::initial_integral(argc, argv);
    if (integral.nonnull())
        sc::Integral::set_default_integral(integral);
    integral = sc::Integral::get_default_integral()->clone();
    return integral;
}

#endif /* DREW_INPUT_H_ */
